{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DoWhy example on the Lalonde dataset\n",
    "\n",
    "Thanks to [@mizuy](https://github.com/mizuy) for providing this example. Here we use the Lalonde dataset and apply IPW estimator to it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: MASS\n",
      "\n",
      "R[write to console]: ## \n",
      "##  Matching (Version 4.9-3, Build Date: 2018-05-03)\n",
      "##  See http://sekhon.berkeley.edu/matching for additional documentation.\n",
      "##  Please cite software as:\n",
      "##   Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching\n",
      "##   Software with Automated Balance Optimization: The Matching package for R.''\n",
      "##   Journal of Statistical Software, 42(7): 1-52. \n",
      "##\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array(['Matching', 'MASS', 'tools', 'stats', 'graphics', 'grDevices',\n",
       "       'utils', 'datasets', 'methods', 'base'], dtype='<U9')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os, sys\n",
    "sys.path.append(os.path.abspath(\"../../../\"))\n",
    "\n",
    "import dowhy\n",
    "from dowhy import CausalModel\n",
    "from rpy2.robjects import r as R\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "#%R install.packages(\"Matching\")\n",
    "%R library(Matching)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%R data(lalonde)\n",
    "%R -o lalonde\n",
    "lalonde = lalonde.astype({'treat':'bool'}, copy=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run DoWhy analysis: model, identify, estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:dowhy.causal_model:Causal Graph not provided. DoWhy will construct a graph based on data inputs.\n",
      "INFO:dowhy.causal_graph:If this is observed data (not from a randomized experiment), there might always be missing confounders. Adding a node named \"Unobserved Confounders\" to reflect this.\n",
      "INFO:dowhy.causal_model:Model to find the causal effect of treatment ['treat'] on outcome ['re78']\n",
      "INFO:dowhy.causal_identifier:Common causes of treatment and outcome:['hisp', 'age', 'nodegr', 'U', 'black', 'educ', 'married']\n",
      "WARNING:dowhy.causal_identifier:If this is observed data (not from a randomized experiment), there might always be missing confounders. Causal effect cannot be identified perfectly.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARN: Do you want to continue by ignoring any unobserved confounders? (use proceed_when_unidentifiable=True to disable this prompt) [y/n] y\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:dowhy.causal_identifier:Instrumental variables for treatment and outcome:[]\n",
      "INFO:dowhy.causal_estimator:INFO: Using Propensity Score Weighting Estimator\n",
      "INFO:dowhy.causal_estimator:b: re78~treat+hisp+age+nodegr+black+educ+married\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Causal Estimate is 1614.1688365342634\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/amit/python-virtual-envs/env3.6/lib/python3.6/site-packages/sklearn/utils/validation.py:73: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  return f(**kwargs)\n"
     ]
    }
   ],
   "source": [
    "model=CausalModel(\n",
    "        data = lalonde,\n",
    "        treatment='treat',\n",
    "        outcome='re78',\n",
    "        common_causes='nodegr+black+hisp+age+educ+married'.split('+'))\n",
    "identified_estimand = model.identify_effect()\n",
    "estimate = model.estimate_effect(identified_estimand,\n",
    "        method_name=\"backdoor.propensity_score_weighting\",\n",
    "        method_params={\"weighting_scheme\":\"ips_weight\"})\n",
    "#print(estimate)\n",
    "print(\"Causal Estimate is \" + str(estimate.value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mnt/c/Users/amit_/code/dowhy/dowhy/interpreters/confounder_distribution_interpreter.py:73: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df[\"weight\"] = df[treated] * (propensity) ** (-1) + (1 - df[treated]) * (1 - propensity) ** (-1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAHwCAYAAABEyLzJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5hcZZX3/e8iAYIQOUbemBCDGjXgg0FCQIxMOAwgHkAFBN/BgDgRjDqKqMg4Eh18ZByVV/BRZIQhzIiCIC+oeOAoIAIGjQgEhigBEgLEAAGUU8h6/ti7odJUd7qTqq7uu7+f6+qrq/beddeqXVWrf7X7rqrITCRJkqRSrdfpAiRJkqR2MvBKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKNiQCb0ScHhH/0qKxJkTEExExoj5/dUR8sBVj1+P9LCJmtmq8flzvSRHxl4h4YKCvu5mIeEtE3LmWl50REYt7WZ8R8eq1HPvNEXFX/Rg4cG3GGAz6+rhdl/thbazLfdMOEXFbRMzo47aLImLvNpdUNHt1n6532PTqNVw2IuI/I+KRiLhpbcboJHtw3wymHtzxwFvfwCcj4vGIeDQiro+IoyPi+doy8+jM/Nc+jtXrzsrMezNzk8x8rgW1z4mI/+42/lszc+66jt3POiYAnwS2y8z/ZyCvuyeZeW1mvrbTdTTxReCb9WPg/+90Me3Wzvuh1QGkHTJz+8y8el3HWZc/7KWwV687e/VqpgN/D4zPzGnN7qMS2IMHTw/ueOCtvSMzRwOvAE4GPgOc2eoriYiRrR5zkJgALM/Mhwb6ipvt00G+n18B3LY2Fxzkt0saCPbqdWOvfsErgEWZ+ddWDFbwY0atkpkd/QEWAXt3WzYNWAW8vj5/NnBSfXor4CfAo8DDwLVUwf2/6ss8CTwBfBqYCCRwFHAvcE3DspH1eFcDXwZuAh4DLga2qNfNABY3qxfYD3gGeLa+vj80jPfB+vR6wOeAe4CHgHOATet1XXXMrGv7C/DPveynTevLL6vH+1w9/t71bV5V13F2k8vOABbX++QhYClwILA/8D/1fjyh2/7/Tb2PlwLfBDZoWJ/AbOAu4O6G8T8DPFDfF6vtO+DlwIV1/XcDH2tYt1F9Hz8C3A58qvt+73Z7EvgY8Od6v/07sF7D+g8AC+rxfgG8ol7+p26PkQ3rui6p98FC4B8bxpkDXAD8d/3Y+GB9P5xZ75clwEnAiB7q7Mt+bHo7gCOAX9eXWQHcAezVcNmr63o2qGv/Xw3rXgb8DRjT5H5YBBwH3FKPex4wqmH9p+ta76/HT+DVTW7bl4DngKfqffnNhtt0dP3YeBT4P0Cs6b5pMv5c4JP16XH1uLPr86+qb3PXvno7ML++vuuBHZr1F6rH2dz6uhfUt3WN+wbYmNWfY09QPW6mAfPqx8aDwNc73U/b+YO92l7d/179DeC++v66GXhLvfwoqt7xXL0vbuzhPuqx3/JCjzwFWE79uGvy+LQH24Or8QdjE62X3wsc06SJfhk4HVi//nlL153ZfSxeaFTn1DtsI5o30SXA6+ttLgT+e01NtD49p2vb7k+ChgfWQuCVwCbAj4D/6lbbf9R1vQF4Gpjcw346h6rBj64v+z/AUT3V2aSJrgQ+X++zf6RqZufW421fP5i2rbffCdgVGFlf1wLg492axGXAFnXtXeP/G1WI7Fq2uN5+Papm93mq5vBKqgazb73+ZKo/hlsA2wC3ruH2JHBVvf2Eel907fMD6n0+ua7/c8D1PT3eqP6wfovqSTWl3i97Nty/z1L9wVmvvl0XAd+heqy8jOqP74d6qLMv+7Gn23FEvU8/Ud9n76VqAF1/4K9u2PZbwL81jPtPwI+bPTbq238TVbPYoq7p6HrdflR/BLcHXkIV9Js22+41dLtNPwE2q2/TMmC/vtw33cb5QMNteB/Vi5XzGtZdXJ/ekSoY7AKMoAoli4ANmzxfTwZ+BWwOjKdqqn3dN6vtx3rZb4DD69ObALt2up+28wd7tb26/736H4At6/o+SdVfRtXrjgCua9i22X3UY7/lhR750Xr8jZpcvz3YHvzC9p1qns2aUrflN1C/imb1JvpFqmbS7BXPamPxQqN6ZZNljU305Ib121G90hzRww5uvPPm0HsTvQL4cMO611IFqK4nX1LNX+pafxNwaJPbNaKuabuGZR8Cru7pgdDt8jOommTXK+PR9XXv0rDNzcCBPVz+48BF3Z5Qe3Yb/xlWf5X6fE1UT4R7u435WeA/69N/pn5C1udnreH2ZLftPwxcUZ/+GfUfl/r8elSvtF/R5P7bhuoV8uiG7b9MfeSlvn+vaVi3NdUfuo0alh0GXNXHx3qz/djT7TiC6hV+4yvzm3jhyd34ONuFKnR0hYl5wCHNHhv17f+HhvNfAU6vT58FfLlh3atZu2Y7veH8+cDxfblvuo3zKqqjAOtRhaYPNTye5gLH1qe/Dfxrt8veCfxdk/v7+T/c9fkP9mPfrLYf62XXAF8AturL/T/Uf7BX26v72aub1PcI8Ib69BH0EnhZQ7+tL39vX6+7l/1jDx4mPXiwzOFtZhzVIfPu/p3qFcovI+LPEXF8H8a6rx/r76F6NbdVn6rs3cvr8RrHHkn1RO7S+E7dv1G9Suluq7qm7mON60cty/OFN388Wf9+sGH9k13XHRGviYifRMQDEfEY8L958f7ovk+XZeZTPVz3K4CX1290eTQiHgVO4IX98HJefB+sSfftX95wXd9ouJ6HgaD5vno58HBmPt5trMZtG6/nFVT3w9KG8b9DdeThRdZiPzbeDoAlWT+re1gPQGbeSPXYmRERr6Nqkpc0q6nW02Ou+/2wpudNf8fv832TmX8C/kp11P0tVEcs7o+I1wJ/R3WUoGvMT3Z7bG1Dk/3Ux9vXl+djl6OA1wB3RMRvI+LtvWxbMnv1C+zVDSLiuIhYEBEr6rE2bVJfT/rSb3t9vNiDXzT+sO7BgzLwRsTOVHfAdd3XZebjmfnJzHwl8E7g2IjYq2t1D0P2tLzLNg2nJ1C9sv8L1Z39koa6RlDNyenruPdTPRgax17J6s2rL/5S19R9rCX9HKevvk01X2lSZr6UquFFt2263/be9sV9wN2ZuVnDz+jM3L9ev5QX3wdr0n37+xuu60PdrmujzLy+yRj3A1tExOhuYzXu18bbdR/VEYetGsZ+aWZu30ONfdmPPd0OgHEREb2sbzSX6t+HhwMX9PIHrTdLqf7N1Ky2Ztb0+O+uP/cNVA31IKo5d0vq8zOp/h02v2HML3Ub8yWZ+f0m4/X39jV60W3NzLsy8zCqP8D/BlwQERv3Y8whz179IvbqWkS8hWqO5iHA5pm5GdWUgO719VRXX/rtmu5Xe/DqhnUPHlSBNyJeWif0H1D9a+OPTbZ5e0S8un4QrqD6l/SqevWDVHOO+usfImK7iHgJ1b/hLqhfYf8PMCoi3hYR61PNd9mw4XIPAhMbP5anm+8Dn4iIbSNiE6pXl+dl5sr+FFfXcj7wpYgYHRGvAI6lmt/TDqOpJoE/Ub9aPWYdx7sJeDwiPhMRG0XEiIh4ff3HEqrb9tmI2DwixlPNyVqTT9Xbb0M1X+q8evnp9VjbA0TEphFxcLMBMvM+qgn2X46IURGxA9Urxqb7NTOXAr8EvlY/VteLiFdFxN/1UGNf9mNPtwOqJ/HHImL9+jZMBi7t4br+G3gXVcM9p4dt1uR84MiImFw/F9b0ear9fb71+b6p/Qr4CNW/raD6991HqP4N2nUE7D+AoyNil6hsXD9fR794uNUeZ+PqsfrqQWDLiNi0a0FE/ENEjMnMVVRv1oAXelHR7NXN2atfVNtKqjmkIyPi88BLe9l+tftoLfptTzXYg18wrHvwYAm8P46Ix6leKfwz8HXgyB62nQRcTvUuvd8A38rMq+p1XwY+F9Vh9eP6cf3/RTX37AGqNy99DCAzV1DN6fku1Sv0v1K9w7XLD+vfyyPid03GPase+xqqd7s+Rd/CXDMfra//z1RHU86tx2+H46gmqT9O9WA+r/fNe1c/Md5O9a+Ru6mOgnyX6t9bUM3Buade90uqfbYmF1PNZZsP/JT6o5Ey8yKqV3o/iOpfWLcCb+1lnMOo5ujdT/UGiRMz8/Jetn8/1Zs5bqea33QBMLaHbfuyH5vejtqNVI/3v1C9I/egzFze7Irq8P47qlfB1/ZSf48y82fAqVRv4lhINTcTqqMszXwDOCiqD44/tQ/j9/e++RXVH6yuZnsd1VG8rvNk5jyqN/Z8k+r+WEg1966ZL1I9f++m6iEX9HLbutd+B1Uo+nPdX15O9QaT2yLiCap9cWhmPtnbOAWwV6+ZvbryC+DnVC9G7qHap739i77ZfdSfftuMPXj18Yd1D+6aYC1pgEVEUv2rbWGTdUdQvRlhej/GOwu4PzM/16L6JlM1xA37e6RrKIiIY6gaZH+OGEkqhD24swa6Bw+WI7yS1kFETATezTp+CUBEvCsiNoyIzamOBPy4lEYbEWOj+mrp9aJ648UnqY7qS9I6sQevWad7sIFXGuIi4l+pjgL8e2bevY7DfYjqMxX/RDXncl3nBA4mG1C9y/tx4Eqqf2V+q6MVSRry7MF91tEe7JQGSZIkFc0jvJIkSSrayE4XsC622mqrnDhxYqfLkFSQm2+++S+ZOWbNW6o7e7KkdmhFXx7SgXfixInMmzev02VIKkhE9OWb/tSEPVlSO7SiLzulQZIkSUUz8EqSJKloBl5JkiQVbUjP4ZWGi2effZbFixfz1FNPdbqUYowaNYrx48ez/vrrd7oUSUOMPbk92tmXDbzSELB48WJGjx7NxIkTiYhOlzPkZSbLly9n8eLFbLvttp0uR9IQY09uvXb3Zac0SEPAU089xZZbbmljbZGIYMstt/TojKS1Yk9uvXb3ZQOvNETYWFvL/SlpXdhDWq+d+9TAK0mSpKI5h1cagiYe/9OWjrfo5LeteZtFi3j729/Orbfe+vyyOXPmsMkmm3Dcccc1vcz8+fO5//772X///VtS59qON2PGDL761a8yderUltQhSY3syYO/J3uEV1LbzJ8/n0svvbTpupUrV7Z0PElS74ZzTzbwSlpnM2bM4DOf+QzTpk3jNa95Dddeey3PPPMMn//85znvvPOYMmUK5513HnPmzOHwww/nzW9+M4cffjjLli3jPe95DzvvvDM777wzv/71rwG46aabeNOb3sSOO+7Ibrvtxp133tl0vL/+9a984AMfYNq0aey4445cfPHFADz55JMceuihTJ48mXe96108+eSTndw9kjSg7Mkv5pQGSS2xcuVKbrrpJi699FK+8IUvcPnll/PFL36RefPm8c1vfhOo/t12++23c91117HRRhvxvve9j0984hNMnz6de++9l3333ZcFCxbwute9jmuvvZaRI0dy+eWXc8IJJ3DhhRe+aLwTTjiBPffck7POOotHH32UadOmsffee/Od73yHl7zkJSxYsIBbbrmFN77xjZ3cNZI04OzJqzPwSuqTnt4927X83e9+NwA77bQTixYt6nGcd77znWy00UYAXH755dx+++3Pr3vsscd44oknWLFiBTNnzuSuu+4iInj22WebjvXLX/6SSy65hK9+9atA9VFB9957L9dccw0f+9jHANhhhx3YYYcd+ndjJWmQsyf3j4FXUp9sueWWPPLII6ste/jhh5//gPANN9wQgBEjRvQ6F2zjjTd+/vSqVau44YYbGDVq1GrbfOQjH2GPPfbgoosuYtGiRcyYMaPpWJnJhRdeyGtf+9q1uUmSNGTZk/vHObyS+mSTTTZh7NixXHnllUDVWH/+858zffr0Hi8zevRoHn/88R7X77PPPpx22mnPn58/fz4AK1asYNy4cQCcffbZPY637777ctppp5GZAPz+978HYPfdd+fcc88F4NZbb+WWW27pz02VpEHPntw/HuGVhqC+fGRNO5xzzjnMnj2bY489FoATTzyRV73qVT1uv8cee3DyySczZcoUPvvZz75o/amnnsrs2bPZYYcdWLlyJbvvvjunn346n/70p5k5cyYnnXQSb3vb23oc71/+5V/4+Mc/zg477MCqVavYdttt+clPfsIxxxzDkUceyeTJk5k8eTI77bRT63eGJNXsyYO/J0dXCh+Kpk6dmvPmzet0GVLbLViwgMmTJ3e6jOI0268RcXNm+oG9a8GerOHCntw+7erLTmmQJElS0Qy8kiRJKtqwm8Pb6q//67ROzRuSpFawJ0saCB7hlSRJUtEMvJIkSSrasJvSIElS28zZtNMVtNacFZ2uQGoJA680FLX6j+oa/qgtX76cvfbaC4AHHniAESNGMGbMGABuuukmNthgg35f5dVXX80GG2zAbrvt1q/LTZw4kXnz5rHVVlv1+zolqS3syYO+Jxt4Ja3Rlltu+fw37syZM4dNNtmE44477vn1K1euZOTI/rWTq6++mk022aTfzVWShjt7cv85h1fSWjniiCM4+uij2WWXXfj0pz/Nn/70J/bbbz922mkn3vKWt3DHHXcA8OMf/5hddtmFHXfckb333psHH3yQRYsWcfrpp3PKKacwZcoUrr32WpYtW8Z73vMedt55Z3beeWd+/etfA9WRjH322Yftt9+eD37wgwzlL8uRpHaxJ/fOI7yS1trixYu5/vrrGTFiBHvttRenn346kyZN4sYbb+TDH/4wV155JdOnT+eGG24gIvjud7/LV77yFb72ta9x9NFHr3ZU4n3vex+f+MQnmD59Ovfeey/77rsvCxYs4Atf+ALTp0/n85//PD/96U8588wzO3yrJWlwsif3zMAraa0dfPDBjBgxgieeeILrr7+egw8++Pl1Tz/9NFA14Pe+970sXbqUZ555hm233bbpWJdffjm333778+cfe+wxnnjiCa655hp+9KMfAfC2t72NzTffvI23SJKGLntyzwy8ktbaxhtvDMCqVavYbLPNnp9T1uijH/0oxx57LO985zu5+uqrmTNnTtOxVq1axQ033MCoUaPaWbIkFcue3DPn8EpaZy996UvZdttt+eEPfwhAZvKHP/wBgBUrVjBu3DgA5s6d+/xlRo8ezeOPP/78+X322YfTTjvt+fNdjXr33Xfn3HPPBeBnP/sZjzzySHtvjCQNcfbkF/MIrzQUDcLPxvze977HMcccw0knncSzzz7LoYceyhve8AbmzJnDwQcfzOabb86ee+7J3XffDcA73vEODjroIC6++GJOO+00Tj31VGbPns0OO+zAypUr2X333Tn99NM58cQTOeyww9h+++3ZbbfdmDBhQodvqSR1Y08e9GKovLuumalTp+a8efP6dRm/t11D0YIFC5g8eXKnyyhOs/0aETdn5tQOlTSk2ZNh0aj3dbqE1hqEQW4wsCe3T7v6slMaJEmSVDQDryRJkopm4JWGiKE8/Wgwcn9KWhf2kNZr5z418EpDwKhRo1i+fLkNtkUyk+XLlxfzcTuSBpY9ufXa3Zf9lAZpCBg/fjyLFy9m2bJlnS6lGKNGjWL8+PGdLkPSEGRPbo929mUDrzQErL/++j1+G47UFxExCrgG2JCq91+QmSdGxNnA3wFdb8c/IjPnR0QA3wD2B/5WL//dwFcuDT725KHHwCtJw8PTwJ6Z+URErA9cFxE/q9d9KjMv6Lb9W4FJ9c8uwLfr35I05DiHV5KGgaw8UZ9dv/7pbQLiAcA59eVuADaLiLHtrlOS2sHAK0nDRESMiIj5wEPAZZl5Y73qSxFxS0ScEhEb1svGAfc1XHxxvaz7mLMiYl5EzHM+o6TBysArScNEZj6XmVOA8cC0iHg98FngdcDOwBbAZ/o55hmZOTUzp44ZM6blNUtSKxh4JWmYycxHgauA/TJzaT1t4WngP4Fp9WZLgG0aLja+XiZJQ46BV5KGgYgYExGb1ac3Av4euKNrXm79qQwHArfWF7kEeH9UdgVWZObSDpQuSevMT2mQpOFhLDA3IkZQHew4PzN/EhFXRsQYIID5wNH19pdSfSTZQqqPJTuyAzVLUksYeCVpGMjMW4Admyzfs4ftE5jd7rokaSA4pUGSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKlrbAm9EbBMRV0XE7RFxW0T8U718i4i4LCLuqn9vXi+PiDg1IhZGxC0R8cZ21SZJkqTho51HeFcCn8zM7YBdgdkRsR1wPHBFZk4CrqjPA7wVmFT/zAK+3cbaJEmSNEy0LfBm5tLM/F19+nFgATAOOACYW282FziwPn0AcE5WbgA2i4ix7apPkiRJw8OAzOGNiIlU3+F+I7B1Zi6tVz0AbF2fHgfc13CxxfWy7mPNioh5ETFv2bJlbatZkiRJZWh74I2ITYALgY9n5mON6zIzgezPeJl5RmZOzcypY8aMaWGlkiRJKlFbA29ErE8Vdr+XmT+qFz/YNVWh/v1QvXwJsE3DxcfXyyRJkqS11s5PaQjgTGBBZn69YdUlwMz69Ezg4obl768/rWFXYEXD1AdJkiRprYxs49hvBg4H/hgR8+tlJwAnA+dHxFHAPcAh9bpLgf2BhcDfgCPbWJskSZKGibYF3sy8DogeVu/VZPsEZrerHkmSJA1PftOaJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0nDQESMioibIuIPEXFbRHyhXr5tRNwYEQsj4ryI2KBevmF9fmG9fmIn65ekdWHglaTh4Wlgz8x8AzAF2C8idgX+DTglM18NPAIcVW9/FPBIvfyUejtJGpIMvJI0DGTlifrs+vVPAnsCF9TL5wIH1qcPqM9Tr98rImKAypWkljLwStIwEREjImI+8BBwGfAn4NHMXFlvshgYV58eB9wHUK9fAWzZZMxZETEvIuYtW7as3TdBktaKgVeShonMfC4zpwDjgWnA61ow5hmZOTUzp44ZM2ada5SkdjDwStIwk5mPAlcBbwI2i4iR9arxwJL69BJgG4B6/abA8gEuVZJawsArScNARIyJiM3q0xsBfw8soAq+B9WbzQQurk9fUp+nXn9lZubAVSxJrTNyzZtIkgowFpgbESOoDnacn5k/iYjbgR9ExEnA74Ez6+3PBP4rIhYCDwOHdqJoSWoFA68kDQOZeQuwY5Plf6aaz9t9+VPAwQNQmiS1nVMaJEmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKlobQu8EXFWRDwUEbc2LJsTEUsiYn79s3/Dus9GxMKIuDMi9m1XXZIkSRpe2nmE92xgvybLT8nMKfXPpQARsR1wKLB9fZlvRcSINtYmScNKRGwTEVdFxO0RcVtE/FO93AMRkoo3sl0DZ+Y1ETGxj5sfAPwgM58G7o6IhcA04DdtKk+ShpuVwCcz83cRMRq4OSIuq9edkplfbdy424GIlwOXR8RrMvO5Aa1aklqgbYG3Fx+JiPcD86ia7yPAOOCGhm0W18teJCJmAbMAJkyY0OZSh4A5m3a6gtaas6LTFUhFysylwNL69OMRsYAe+mzNAxGSijHQb1r7NvAqYApV4/1afwfIzDMyc2pmTh0zZkyr65Ok4tX/fdsRuLFe9JGIuKV+78Xm9bJxwH0NF+vxQIQkDXYDGngz88HMfC4zVwH/QXW0AGAJsE3DpuPrZZKkFoqITYALgY9n5mOs44GIiJgVEfMiYt6yZctaXq8ktcKABt6IGNtw9l1A1yc4XAIcGhEbRsS2wCTgpoGsTZJKFxHrU4Xd72Xmj2DdD0T4XzdJQ0Hb5vBGxPeBGcBWEbEYOBGYERFTgAQWAR8CyMzbIuJ84HaqN1bM9o0RktQ6ERHAmcCCzPx6w/Kx9fxeePGBiHMj4utUb1rzQISkIaudn9JwWJPFZ/ay/ZeAL7WrHkka5t4MHA78MSLm18tOAA7zQISk0nXiUxokSQMsM68DosmqS3u5jAciJBXBrxaWJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKL1KfBGxJv7skyS1H72ZEnqn74e4T2tj8skSe1nT5akfhjZ28qIeBOwGzAmIo5tWPVSYEQ7C5MkrW5denJEbAOcA2wNJHBGZn4jIrYAzgMmAouAQzLzkYgI4BvA/sDfgCMy83etvUWSNDDWdIR3A44swq8AABGXSURBVGATqmA8uuHnMeCg9pYmSepmXXrySuCTmbkdsCswOyK2A44HrsjMScAV9XmAtwKT6p9ZwLdbe1MkaeD0eoQ3M38F/Coizs7MewaoJklSE+vSkzNzKbC0Pv14RCwAxgEHADPqzeYCVwOfqZefk5kJ3BARm0XE2HocSRpSeg28DTaMiDOo/uX1/GUyc892FCVJ6tU69eSImAjsCNwIbN0QYh+gmvIAVRi+r+Fii+tlqwXeiJhFdQSYCRMm9O9WSNIA6Wvg/SFwOvBd4Ln2lSNJ6oO17skRsQlwIfDxzHysmqpbycyMiOzPeJl5BnAGwNSpU/t1WUkaKH0NvCsz0/lbkjQ4rFVPjoj1qcLu9zLzR/XiB7umKkTEWOChevkSYJuGi4+vl0nSkNPXjyX7cUR8OCLGRsQWXT9trUyS1JN+9+T6UxfOBBZk5tcbVl0CzKxPzwQublj+/qjsCqxw/q6koaqvR3i7muGnGpYl8MrWliNJ6oO16clvBg4H/hgR8+tlJwAnA+dHxFHAPcAh9bpLqT6SbCHVx5Id2ZrSJWng9SnwZua27S5EktQ3a9OTM/M6IHpYvVeT7ROY3d/rkaTBqE+BNyLe32x5Zp7T2nIkSWtiT5ak/unrlIadG06Pojoa8Duqb+2RJA0se7Ik9UNfpzR8tPF8RGwG/KAtFUmSemVPlqT+6eunNHT3V8B5vZI0ONiTJakXfZ3D+2OqdwADjAAmA+e3qyhJUs/syZLUP32dw/vVhtMrgXsyc3Eb6pEkrZk9WZL6oU9TGjLzV8AdwGhgc+CZdhYlSeqZPVmS+qdPgTciDgFuAg6m+lDyGyPioHYWJklqzp4sSf3T1ykN/wzsnJkPAUTEGOBy4IJ2FSZJ6pE9WZL6oa+f0rBeV2OtLe/HZSVJrWVPlqR+6OsR3p9HxC+A79fn30v1PeuSpIFnT5akfug18EbEq4GtM/NTEfFuYHq96jfA99pdnCTpBfZkSVo7azrC+/8BnwXIzB8BPwKIiP9Vr3tHW6uTJDWyJ0vSWljTnK+tM/OP3RfWyya2pSJJUk/syZK0FtYUeDfrZd1GrSxEkrRG9mRJWgtrCrzzIuIfuy+MiA8CN7enJElSD+zJkrQW1jSH9+PARRHx//JCM50KbAC8q52FSZJexJ4sSWuh18CbmQ8Cu0XEHsDr68U/zcwr216ZJGk19mRJWjt9+hzezLwKuKrNtUiS+sCeLEn94zfzSJIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFa1vgjYizIuKhiLi1YdkWEXFZRNxV/968Xh4RcWpELIyIWyLije2qS5IkScNLO4/wng3s123Z8cAVmTkJuKI+D/BWYFL9Mwv4dhvrkiRJ0jDStsCbmdcAD3dbfAAwtz49FziwYfk5WbkB2CwixrarNkmSJA0fAz2Hd+vMXFqffgDYuj49DrivYbvF9bIXiYhZETEvIuYtW7asfZVKkiSpCB1701pmJpBrcbkzMnNqZk4dM2ZMGyqTJElSSQY68D7YNVWh/v1QvXwJsE3DduPrZZIkSdI6GejAewkwsz49E7i4Yfn7609r2BVY0TD1QZIkSVprI9s1cER8H5gBbBURi4ETgZOB8yPiKOAe4JB680uB/YGFwN+AI9tVlyRJkoaXtgXezDysh1V7Ndk2gdntqkWSJEnDl9+0JkmSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JWkYSAizoqIhyLi1oZlcyJiSUTMr3/2b1j32YhYGBF3RsS+nalaklrDwCtJw8PZwH5Nlp+SmVPqn0sBImI74FBg+/oy34qIEQNWqSS1mIFXkoaBzLwGeLiPmx8A/CAzn87Mu4GFwLS2FSdJbWbglaTh7SMRcUs95WHzetk44L6GbRbXy14kImZFxLyImLds2bJ21ypJa8XAK0nD17eBVwFTgKXA1/o7QGaekZlTM3PqmDFjWl2fJLWEgVeShqnMfDAzn8vMVcB/8MK0hSXANg2bjq+XSdKQZOCVpGEqIsY2nH0X0PUJDpcAh0bEhhGxLTAJuGmg65OkVhnZ6QIkSe0XEd8HZgBbRcRi4ERgRkRMARJYBHwIIDNvi4jzgduBlcDszHyuE3VLUisYeCVpGMjMw5osPrOX7b8EfKl9FUnSwHFKgyRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JWmYiIizIuKhiLi1YdkWEXFZRNxV/968Xh4RcWpELIyIWyLijZ2rXJLWjYFXkoaPs4H9ui07HrgiMycBV9TnAd4KTKp/ZgHfHqAaJanlDLySNExk5jXAw90WHwDMrU/PBQ5sWH5OVm4ANouIsQNTqSS1loFXkoa3rTNzaX36AWDr+vQ44L6G7RbXy1YTEbMiYl5EzFu2bFl7K5WktWTglSQBkJkJZD8vc0ZmTs3MqWPGjGlTZZK0bgy8kjS8Pdg1VaH+/VC9fAmwTcN24+tlkjTkGHglaXi7BJhZn54JXNyw/P31pzXsCqxomPogSUPKyE4XIEkaGBHxfWAGsFVELAZOBE4Gzo+Io4B7gEPqzS8F9gcWAn8DjhzwgiWpRQy8kjRMZOZhPazaq8m2Ccxub0WSNDCc0iBJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKJ15JvWImIR8DjwHLAyM6dGxBbAecBEYBFwSGY+0on6JEmSVI5OHuHdIzOnZObU+vzxwBWZOQm4oj4vSZIkrZPBNKXhAGBufXoucGAHa5EkSVIhOhV4E/hlRNwcEbPqZVtn5tL69APA1p0pTZIkSSXpyBxeYHpmLomIlwGXRcQdjSszMyMim12wDsizACZMmND+SiVJkjSkdeQIb2YuqX8/BFwETAMejIixAPXvh3q47BmZOTUzp44ZM2agSpYkSdIQNeCBNyI2jojRXaeBfYBbgUuAmfVmM4GLB7o2SZIklacTUxq2Bi6KiK7rPzczfx4RvwXOj4ijgHuAQzpQmyRJkgoz4IE3M/8MvKHJ8uXAXgNdjyRJkso2mD6WTJIkSWo5A68kSZKK1qmPJZMkSdKazNm00xW01pwVHblaj/BKkiSpaB7hlSRJxZh4/E87XUJLLRrV6QrK4BFeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKlofvGEtI6K+5Dzk9/W6RIkSWopj/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloftOapNXN2bTTFbTWnBWdrkCS1GEe4ZUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNT2mQpGEuIhYBjwPPASszc2pEbAGcB0wEFgGHZOYjnapRktaFR3glSQB7ZOaUzJxanz8euCIzJwFX1OclaUgy8EqSmjkAmFufngsc2MFaJGmdGHglSQn8MiJujohZ9bKtM3NpffoBYOvOlCZJ6845vJKk6Zm5JCJeBlwWEXc0rszMjIhsdsE6IM8CmDBhQvsrlaS14BFeSRrmMnNJ/fsh4CJgGvBgRIwFqH8/1MNlz8jMqZk5dcyYMQNVsiT1i4FXkoaxiNg4IkZ3nQb2AW4FLgFm1pvNBC7uTIWStO6c0iBJw9vWwEURAdXfhHMz8+cR8Vvg/Ig4CrgHOKSDNUrSOjHwStIwlpl/Bt7QZPlyYK+Br0iSWs8pDZIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEMvJIkSSqagVeSJElFM/BKkiSpaAZeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkohl4JUmSVDQDryRJkopm4JUkSVLRDLySJEkqmoFXkiRJRTPwSpIkqWgGXkmSJBXNwCtJkqSiGXglSZJUNAOvJEmSimbglSRJUtEGXeCNiP0i4s6IWBgRx3e6HkkazuzJkkowqAJvRIwA/g/wVmA74LCI2K6zVUnS8GRPllSKQRV4gWnAwsz8c2Y+A/wAOKDDNUnScGVPllSEkZ0uoJtxwH0N5xcDuzRuEBGzgFn12Sci4s4Bqm1QCtgK+Eun62iZL0SnKxj2fEzxinaUMUTZk/vJ549azccU0IK+PNgC7xpl5hnAGZ2uY7CIiHmZObXTdagcPqbUH/bk1fn8Uav5mGqNwTalYQmwTcP58fUySdLAsydLKsJgC7y/BSZFxLYRsQFwKHBJh2uSpOHKniypCINqSkNmroyIjwC/AEYAZ2XmbR0ua7DzX4lqNR9TAuzJa8nnj1rNx1QLRGZ2ugZJkiSpbQbblAZJkiSppQy8kiRJKpqBd4jy6z7VahFxVkQ8FBG3droWaSiyL6uV7MmtZeAdgvy6T7XJ2cB+nS5CGorsy2qDs7Ent4yBd2jy6z7Vcpl5DfBwp+uQhij7slrKntxaBt6hqdnXfY7rUC2SJPuyNKgZeCVJklQ0A+/Q5Nd9StLgYl+WBjED79Dk131K0uBiX5YGMQPvEJSZK4Gur/tcAJzv131qXUXE94HfAK+NiMURcVSna5KGCvuyWs2e3Fp+tbAkSZKK5hFeSZIkFc3AK0mSpKIZeCVJklQ0A68kSZKKZuCVJElS0Qy8GrYi4uURcUE/L3N2RBzUrpokabiyJ6udDLwaFiJiZPfzmXl/ZtooJWmA2ZM10Ay8GtQiYmJE3FG/iv+fiPheROwdEb+OiLsiYlr985uI+H1EXB8Rr60ve0REXBIRVwJXNDk/MSJurbcdERH/HhG/jYhbIuJD9fKIiG9GxJ0RcTnwso7tDEnqMHuyhqqRa95E6rhXAwcDH6D6+s73AdOBdwInAO8H3pKZKyNib+B/A++pL/tGYIfMfDgijuh2fmLDdRwFrMjMnSNiQ+DXEfFLYEfgtcB2wNbA7cBZbbytkjTY2ZM15Bh4NRTcnZl/BIiI24ArMjMj4o/ARGBTYG5ETAISWL/hspdl5sO9nO+yD7BDw1ywTYFJwO7A9zPzOeD++kiEJA1n9mQNOQZeDQVPN5xe1XB+FdVj+F+BqzLzXfURgqsbtv9rt7G6n+8SwEcz8xerLYzYf+1KlqRi2ZM15DiHVyXYFFhSnz5iLcf4BXBMRKwPEBGviYiNgWuA99bzycYCe6xrsZJUOHuyBh0Dr0rwFeDLEfF71v6/Ft+lmgv2u/pNE9+px7oIuKtedw7wm3UvV5KKZk/WoBOZ2ekaJEmSpLbxCK8kSZKKZuCVJElS0Qy8kiRJKpqBV5IkSUUz8EqSJKloBl5JkiQVzcArSZKkov1fYz/cUbiJQP4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "estimate.interpret(method_name=\"confounder_distribution_interpreter\",var_type='discrete',\n",
    "                   var_name='married', fig_size = (10, 7), font_size = 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sanity check: compare to manual IPW estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Causal Estimate is 1639.8546888423534\n"
     ]
    }
   ],
   "source": [
    "df = model._data\n",
    "ps = df['ps']\n",
    "y = df['re78']\n",
    "z = df['treat']\n",
    "\n",
    "ey1 = z*y/ps / sum(z/ps)\n",
    "ey0 = (1-z)*y/(1-ps) / sum((1-z)/(1-ps))\n",
    "ate = ey1.sum()-ey0.sum()\n",
    "print(\"Causal Estimate is \" + str(ate))\n",
    "\n",
    "# correct -> Causal Estimate is 1634.9868359746906"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
